## Robustness Checks (Results are estimated by lfe package version 2.8-7)
## These results are reported in Table 4-6
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate the short sample (shorter in both pre and post-shutdown)
data4=subset(data3, T>40&T<127)

##############################################################################################################################
########################
## Post shutdown placebo (Table 4)
# Long sample
all_AOD_long=felm(AOD~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data3)
summary(all_AOD_long)
all_SO2_long=felm(SO2_MASS..tons.~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data3)
summary(all_SO2_long)
all_NOX_long=felm(NOX_MASS..tons.~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data3)
summary(all_NOX_long)
# Short sample
all_AOD_short=felm(AOD~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data4)
summary(all_AOD_short)
all_SO2_short=felm(SO2_MASS..tons.~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data4)
summary(all_SO2_short)
all_NOX_short=felm(NOX_MASS..tons.~ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data4)
summary(all_NOX_short)
stargazer(all_SO2_long, all_NOX_long,all_AOD_long,all_SO2_short, all_NOX_short, all_AOD_short)

##############################################################################################################################
##############################################################################################################################
## Time shifting placebo 1 (Table 5)
data4=data3
data5=subset(data4, T>40&T<127)
## Generate placebo treatments
data4$pre=0
data4$pre[which( data4$DATE>="2018-12-01")]=1
data6=subset(data4, T<=68)
all_AOD_pre=felm(AOD~   pre+ ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                 data = data6)
summary(all_AOD_pre)
all_SO2_pre=felm(SO2_MASS..tons.~ pre + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                 data = data6)
summary(all_SO2_pre)
all_NOX_pre=felm(NOX_MASS..tons.~ pre + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                 data = data6)
summary(all_NOX_pre)
all_AOD=felm(AOD~   pre + shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
             data = data4)
summary(all_AOD)
all_SO2=felm(SO2_MASS..tons.~ pre + shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
             data = data4)
summary(all_SO2)
all_NOX=felm(NOX_MASS..tons.~ pre+ shutdown + post + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
             data = data4)
summary(all_NOX)
stargazer(all_SO2_pre, all_NOX_pre, all_AOD_pre, all_SO2, all_NOX, all_AOD)

##############################################################################################################################
##############################################################################################################################
## Time shifting placebo 2 (Table 6)
# Long sample
data4$pre=0
tt="2018-10-22"
for (i in 1:6) {
  tt=as.Date(tt)+7
  data4$pre[which( data4$DATE>=as.Date(tt)-7& data4$DATE<as.Date(tt))]=i
}
tt="2018-12-01"
for (i in 7:10) {
  tt=as.Date(tt)+7
  data4$pre[which( data4$DATE>=as.Date(tt)-7& data4$DATE<as.Date(tt))]=i
}
data4$pre[which( data4$DATE>"2018-12-28")]=0
## set 2018-12-28 as reference
data4$pre[which( data4$DATE=="2018-11-01")]=0
data6=subset(data4, T<=68)
all_AOD_long=felm(AOD~   factor(pre) + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data6)
summary(all_AOD_long)
all_SO2_long=felm(SO2_MASS..tons.~ factor(pre)+ ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data6)
summary(all_SO2_long)
all_NOX_long=felm(NOX_MASS..tons.~ factor(pre)+ ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data6)
summary(all_NOX_long)
# Short sample
data5$pre=0
tt="2018-12-01"
for (i in 1:4) {
  tt=as.Date(tt)+7
  data5$pre[which( data5$DATE>=as.Date(tt)-7& data5$DATE<as.Date(tt))]=i+6
}
data5$pre[which( data5$DATE>"2018-12-28")]=0
## set 2018-12-28 as reference
data5$pre[which( data5$DATE=="2018-12-28")]=0
data6=subset(data5, T<=68)
all_AOD_short=felm(AOD~   factor(pre) + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data6)
summary(all_AOD_short)
all_SO2_short=felm(SO2_MASS..tons.~ factor(pre)+ ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data6)
summary(all_SO2_short)
all_NOX_short=felm(NOX_MASS..tons.~ factor(pre)+ ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data6)
summary(all_NOX_short)
stargazer(all_SO2_short, all_NOX_short, all_AOD_short, all_SO2_long, all_NOX_long, all_AOD_long)

